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Use of High-Resolution 

Upwind Scheme for Vortical Flow Simulations* 
KozoFUJIl** and Shigeru OBAYASHI*** 


ABSTRACT 

For vortical flow simulations at high Reynolds numbers, it is important to keep 
the artificial dissipation as small as possible since it induces unphysical decay of the 
vortex strength. One way to accomplish this is to decrease the grid spacing. Another way 
is to use computational schemes having little dissipation. In the present paper, one of 
the high-resolution upwind schemes called MUSCL with Roe’s average’ is applied to 
vortical flow fields. 

Two examples are considered. One is the leading^dge separation-vortex flow over 
a strake-delta wing. The other is a high-angle of attack supersonic flow over a spaceplane- 
like geometry. Comparison with the central difference solutions indicates that the present 
upwind scheme is less dissipative and thus has better resolution for the vortical flows. 


lC A.l.&'jA' 3 75 > £ <3 A' $ x. 3 • £ i * 3 o ft, 

zzm t> j oift o ',n > hrk-z&h -> e t-j j -ooijmtW'bztjiK 

k'ni wf a- 4 -fijifl-r 5 c £ t* 5 0 i- -> x^ua 

morm'/sttn 

uixu re* ->T t, fj&K a zzk tun 
m cS6 ° c film stir itzfimiA t.mtJi-Roenvmtm'tcMu- 

SCL&£ftJHJ l, t P>C'fi : :5)cDt%&coM£ltfi l ffiji LTIt'/y'ju-r 

* 9 • -o 

CO) F-iicIi, 'Till l 'tVtT'Sr i, ->z i ' 6 tz if A A < 

H ' T b fffrAfttT-’v/j.A c. £ /jiHlj uj ^ £ A" /- o 


1. INTRODUCTION 

The flow over aircraft and missiles at moderate 
to high angles of attack is characterized by the 
presence of large spiral vortices on the leeward 
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side ol the body. These separation vortices in- 
duce low pressure on the body upper surface, 
and this low pressure is the predominant factor 
of the resulting aerodynamic characteristic of 
the body. Research on such flow is of great 
importance practically as well as physically be- 
cause understanding of the separated and vortical 
flow fields may lead to the control of vortex 
behavior and eventually to the enhancement of 
flight vehicle performance. 
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The rapid progress of supercomputers and 
numerical methods has made computer simula- 
tion of such vortical flows feasible. (Recent 
efforts in computational methods and applica- 
tion results were surveyed by Newsome and 
Kandil, see Ref. 1). Three-dimensional Reynolds- 
averaged Navier-Stokes equations used for 
various flow simulations have become a primary 
method for the vortical flow simulations since 
they can describe separated vortical flows with 
no special treatment. There exists one important 
fact that should be kept in mind. Typical com- 
putational grids used for Navier-Stokes simula- 
tions are fine only near the body surface to 
resolve viscous layers. These grid distributions 
are adequate for those flows where important 
phenomena only occur near the body surface as 
embedded shock waves or flow separation. For 
vortical flow simulations, on the other hand, not 
only the region near the body surface but also 
regions away from the body surface are impor- 
tant. Since the strength of vortices which are 
located away from the body is a dominant factor 
of the flow field, grid resolution away from the 
body is critical for an accurate simulation of 
vortical flow. Rai in Ref. 2 pointed out that 
conventional, spatially second-order accurate 
finite-difference schemes are much too dissipative 
for calculations involving vortices that travel 
large distances. Reference 3 studied the effect 
of grid resolution for vortical flow simulations 
and found the important result that vortex 
breakdown phenomenon cannot be predicted 
unless a sufficiently fine grid is used. 

Now, how can we obtain accurate solutions 
for the vortical flows with computer memory 
and time constraints? High-order upwind differ- 
encing has become popular in recently-developed 
Euler methods for compressible inviscid flows. 
This feature lias been extended in the straight- 
forward manner for the evaluation of the con- 
vective terms of Navier-Stokes computations (see 
Refs. 4-6, for instance). Discontinuities are more 
sharply captured by these upwind methods com- 


pared to those by conventional central difference 
method with additional numerical dissipation, 
since high-order upwind methods introduce 
minimum amount of dissipation needed to 
prevent oscillations. Recently, the matrix form 
of the dissipation terms implicitly introduced by 
upwind methods was studied 7 ^ 9 ^ and it was 
shown that such terms in the upwind schemes 
such as Roe’s flux difference splitting become 
small in the viscous layers. As artificial dissipa- 
tion should be kept to a minimum and viscous 
effects near the body surface should be correctly 
evaluated in vortical flow simulations, high- 
resolution upwind method may be adequate for 
vortical flow simulations. 

The object of the present work is to demon- 
strate the capability of high-resolution upwind 
method for accurate vortical flow simulations. 
Two flow fields are considered as application 
examples. One is a subsonic flow over a strake- 
delta wing. Simulation of the same flow fields 
was already conducted by the first author using 
a conventional central difference method. The 
same grid distributions: fine grid and medium 
grid are used for comparison. The second ex- 
ample is a supersonic flow over a spaceplane-like 
geometry. The flow field is much more com- 
plicated in this example since there exist bow 
shock wave, wing shock wave, fuselage vortex, 
strake vortex and wing vortex. 

2. GOVERNING EQUATIONS AND 
NUMERICAL ALGORITHM 

Compressible Navier-Stokes Equations 

The basic equations under consideration are 
the unsteady Navier-Stokes equations written for 
a body-fitted coordinate system (£, t?, f). 

3rC + M +8 ^ + M = Re ~ l M 0) 

In Eq. (1), the thin-layer approximation has 
been introduced. The use of thin-layer Navier- 
Stokes equations is justified because the viscous 
effects are confined to a thin layer near the wall 
and are dominated by the viscous terms asso- 
ciated with the strain rates normal to the wall, 
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and because the flow away from the body is 
essentially rotational inviscid. The contribution 
of the viscous terms in the shear layer rolling up 
from the surface wall and core of the vortices is 
assumed to be negligible. It should be noted that 
the viscous terms are not properly evaluated in 
these regions even by the full Navier-Stokes 
equations because of the g r id deficiency. 

The pressure, density, and velocity compo- 
nents are related to the energy for an ideal gas by 

p = (7- 1) [e -p (n 2 +^ 2 +co 2 )/2] (2) 

LU-ADI Algorithm with Upwind Feature 

The time-integration algorithm used here is 
the LU-ADI factorization method proposed by 
Obayashi et al. for two-dimensional problems. 
An extension to three-dimensional problems is 
described in Ref. 11. This algorithm decomposes 
the usual block tridiagonal system of Beam and 
Warming’s into the product of lower and upper 
scalar bidiagonal matrices using a diagonal form 
and an approximate LDU decomposition. The 
basic algorithm is first order accurate in time. In 
this original LU-ADI scheme, the convective 
terms are evaluated using second- or fourth -order 
central differencing and the viscous terms are 
evaluated using second-order centra] differencing 
in the right-hand side. Since the delta form is 
used, steady state solutions are indifferent to the 
left-hand-side operators, and depend only on the 
right-hand-side steady state description. Thus, 
steady state solutions can be improved by 
modifying the right-hand-side discretization 
method. 

In the right-hand side, convective terms are 
now evaluated using flux difference splitting by 
Roe. 12) The MUSCL interpolation is used for the 
higher-order extension following Ref. 9. Central 
differencing that is adopted in the original LU- 
ADI code is used just for the comparison, In that 
case, nonlinear artificial dissipation terms are 
added (see Ref. 13), and their coefficients are set 
to be reasonable values that have been successful- 
ly used for many practical computations. 


Higher-order extension of flux difference 
splitting using the MUSCL approach is found in 
Ref. 9, but is briefly described again. When the 
convective terms are differenced with the flux- 
difference splitting of Roe, the spatial derivatives 
are written in the conservative form as a flux 
balance. For instance in the ^-direction, 

(ff) = ^i + \l2 ~ ^)-l/2)/(£/+l/2 “ fy-1/2) 

(3) 

The numerical flux F/+ 1/2 can be written as the 
solution to an approximate Rieman problem and 
the necessary metric terms are evaluated at the 
cell interface j + 1/2. 

^>1/2=7 l f(Q l ) + f{Q r ) 

-\A\(Q r ~Q l )] M /2, (4) 

where F is the flux vector and A is the corre- 
sponding Jacobian matrix computed using the 
Roe’s average state Q L and Q R are the state 
variables to the left and right of the half-cell 
interface. These state variables are determined 
from the locally one-dimensional non-oscillatory 
interpolations called MUSCL approach. Primi- 
tive variables q [p, u, v, w, P] are used for that 
purpose, and high-order accurate monotone 
differencing is given by a one-parameter k. 

= £ 7/ + 7 K 1 - ks)A_+(1 +KS)AJ, 

(Qr) !+ \I2 = <7/+i - j [0 - «*) A. + (1 + *s) A J /+ i . 

( 5 ) 

where 

(A + )/ = «/♦ 1 ?/• = q i ' q i \ • 

2A + A_ + e 
* = (AJ 5 + (A_) 2 + e 

s is the Van Alabama limiter and k is a small 
constant to prevent zero division. For all the 
results here, third-order accuracy corresponding 
to k = 1/3 is used. Near the boundary, the 
MUSCL interpolation goes down to the first- 
order. 
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3. RESULTS 

Computed Results for Subsonic Flow over 
a Strake-Delta Wing 

First example is the subsonic flow over a 
strake-delta wing. The flow field was extensively 
studied experimentally by Brennenstuhl and 
Hummel 14 -* in a low-speed wind tunnel and com- 
putationally by one of the present authors. The 
freestream Mach number is 0.3, and the Reynolds 
number based on the root chord is 1.3 x 10 6 in 
the following computations. 

The angle of attack for the first case is 12 
degrees. At this angle of attack, there exist two 
vortices over the upper surface of the wing; one 
emanating from the strake leading edge and the 
other from the main-wing leading edge. These 
two vortices merge together over the main-wing 
surface because of the mutual interaction. 
Figures la and lb show the overall view of the 
spanwise total pressure contour plots at several 
chordwise stations. The result obtained by the 
upwind differencing is plotted in Fig. la, and 
the result by the central (with added dissipation) 
differencing is plotted in Fig. lb. The total 
number of grid points is about 850,000; 119 
points in the chordwise direction, 101 points 
circumferentially and 71 points in the normal 
direction. Details of the grid generation and the 
grid distribution can be found in Ref. 15. Both 
results indicate the existence of two vortices over 
the wing surface and their interaction. It seems 
that the merging of the two vortices is slightly 
delayed in the upwind solution. The correspond- 
ing particle path traces showing the vortex 
trajectories are shown in Fig. 2 for the upwind 
result. The comparison of the computed vortex 
trajectories with the experiment that was already 
done for the central differencing result in Ref. 15 
is also presented in this figure. It is clear that 
merging of two vortices are delayed in the up- 
wind result, hut still earlier than the experiment 
at the same Reynolds number. For further com- 
parison. the spanwise total pressure contour plots 


(with twice the number of contours lines) and 
density contour plots as 65% chordwise location 
are depicted in Figs. 3a-4b. Figures 3a and 3b are 
the results obtained by the upwind differencing, 
and Figs. 4a and 4b are the results by the central 
differencing. Both contour plots show strong 
gradients due to the viscous layers near the sur- 
face. Upwind result shows that this region is 
thinner compared to the central differencing 
result. This may be reasonable since numerical 
dissipation introduced by the present upwindlng 
is small in the viscous boundary layers as ex- 
plained in Ref. 9. 

The same computation was carried out using 
smaller number of grid points (called medium 
grid in Ref. 3, about 120,000 in total). Com- 
pared to the previous grid, the number of the 
grid points are decreased in all the directions. 
Figures 5a and 5b represent the total pressure 
contour plots obtained by the upwind and 
central difference computations, respectively. 
The upwind result shown in Fig. 5a indicates the 
existence of two vortices and their merging 
process although the inner vortex is not as 
distinct as the fine grid result. On the other 
hand, the central difference result in Fig. 5b 
shows only one flattened vortex instead ot two 
vortices. 

It is recognized from these results that the 
present upwind scheme has better resolution 
than the conventional central difference scheme 
on the same grid although grid resolution itself is, 
of course, an important tactor lor an accurate 
How simulation. Upwind scheme is more "vortex- 
preserving" than central diflerencing scheme 
(with added dissipation) since it has a lower level 
of dissipation. In the present upwind scheme 
where the dissipation terms are constructed in 
the matrix form, each characteristic wave has its 
own minimum dissipation. On the other hand, 
central difference scheme where dissipation 
terms are constructed in the scalar lorm, requires 
amount of dissipation which is large enough lor 
all the waves. This is the reason that the upwind 
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a) upwind difference result 



b) central difference result 


Fig, 1 Spanwise total pressure contour plots; a - 12 
overall view of the fine grid solution - 
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a) total pressure contour 



b) density contour 


Fig. 3 Spanwise contour plots at 65 % chord by 
the upwind scheme: a = 1 2° 

scheme like the one used here is less dissipative 
and has better resolution. 

Of course, solution of both central and 
upwind difference schemes should converge to 
the solution of the original partial differential 
equations as the computational grid is refined. 
The accuracy estimation based on the idea ot 
the Tailor expansion is important but not good 
enough for the system of nonlinear equations. 
What we need in numerical schemes is the better 
representation of the properties of original 
partial differential equations and, in that sense, 
upwind difference scheme shows better result 
than that of the central difference scheme for 
the grid distributions feasible under the memory 
restriction of the current supercomputers. 

The angle of attack for the second case is 30 
degrees. Only the calculations for the medium 
grid (previously mentioned grid ot about 1 20,000 
points) computation is carried out. At this angle 



a) total pressure contour 



Fig. 4 Spanwise contour plots at 65% chord by 
the central difference scheme: a = 12° 

of attack, both the experiment and the com- 
putational results using the fine grid showed that 
vortex breakdown takes place near the trailing 
edge. The computed total pressure contour plots 
are presented in Figs. 6a and 6b. An abrupt in- 
crease of the vortex-core is observed near the 
trailing edge in the upwind result shown in Fig. 
5a. This indicates that the vortex has undergone 
breakdown. In fact, the plot ot the streamwise 
velocity (although not shown here) showed that 
there exists the reverse flow region near the 
trailing edge. The central difference result shown 
in Fig. 6b, on the other hand, does not show 
such a sudden change. Again, the resolution is 
enhanced by the use of the present upwind 
scheme at least on the grid used here although a 
slight increase of the number of grid points may 
introduce breakdown phenomenon also in the 
central difference result. 
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a) upwind result 



hi central difference result 


Fig. 6 Computed total pressure contour plots: a = 30 
overall view of 1 the medium grid solution 


Computed Results for Supersonic Flow over 
a Spaceplane 

Computations were carried out at the Mach 
number, 1.5, the Reynolds number, 1.0 x 10* 
based on the maximum span length, and the 
angle of attack 20 degrees. The configuration has 
only fuselage and strake-wing corresponding to 
the experimental model. Figure 7 shows an over- 
view o! the surface grid over this geometry. 
Topologically, the grid is of C-type in the chord- 


wise direction and 0-type in the circumferential 
direction, and consists of 89 points in the Jtord- 
wise direction, 103 points in the circumferential 
direction and 40 points from the body to the 
outer boundary. The minimum spacing near 
the body surface is carefully set to be 0 (10' s ) 
with the span length unity. The outer boundary 
is located outside of the bow shock which is 
captured. In the computation, bilateral sym- 
metry is assumed and the half of the volume is 


i* 
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Fig. 7 Overall view of the computational grid for the spaeeplane 


solved. At the end of the body geometry, zero- 
th order extrapolation is applied as the outflow 
condition. Even for this simplified geometry, 
the flow field is complicated since there exist 
bow shock, strake and wing shocks and leading- 
edge separation vortices. Thus, it is important 
to use an accurate solution scheme both for 
linear waves and nonlinear waves. Figures 8 a 
and 8b show the overall density contour plots 
for the upwind and central difference solutions, 
respectively. Since the present upwind scheme 
is TV D- type, all the shock waves; bow shock 
wave from the nose, wing shocks from the strake 
and the wing, are captured better in Fig. 8a com- 
pared to the central difference result shown in 
Fig. 8b. Stream wise periodic structure observed 
near the symmetry' plane on the fuselage surface 
indicates the separation vortex from the fuselage 
is developed in the stream wise direction. It is 
noticed that contours over the surface especially 
near the wing-fuselage junction are significantly 
different between Figs 8a and Hb. Also noticed 
is the kink of the contour lines on the fuselage 
surface where the strake begins. This may he due 
to file distribution of the grid and the effect is 
pronounced when the upwind differencing is 
used. 

Two cross-sectional plots are presented to 
show the detail. C rossflow density plots at the 
forward fuselage section j- 1 W7< chord) are 
shown in Fig. l >, with the upwind result presented 
“ii the left, arid the central difference result on 
the right. Again, bow shock is clearer in the up- 


wind result. The central difference solution 
demonstrates stronger primary' vortex although 
both solutions produce the primary vortex. This 
may be caused by the artificial dissipation terms 
implemented in the central difference computa- 
tion. In other words, central difference solution 
represents sort of relatively “lower Reynolds 
number flow” near the body surface because of 
the added dissipation and thus the flow tends 
to separate more easily. The cross-sectional 
geometry is composed of the flat region in the 
bottom and circular region on the top. At the 
junction, another separation may lx* recognized 
in the upwind solution although not clearly seen. 
Figure 10 shows the similar plots for the rear 
portion of the geometry f % chord). One 

interesting feature here is the existence of the 
crossflow shock wave in the upwind solution. 
Experimental tlow visualization for the crossflow 
was not carried out and thus it is impossible to 
say which solution is physically correct. How- 
ever, it is important that different flow fields are 
obtained h\ the two different solution schemes. 

4 CONCLUSION 

For vortical flow simulations at high Reynolds 
number, it is important to keep artificial dissipa- 
tion as small as possible since resolution is critical 
for an accurate simulation. One way is to de- 
crease the truncation error by reducing the grid 
spacing. Numerical experiments shown in the 
present paper indicated that the present upwind 
scheme Fiu> better resolution on the same grid 
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a ) upwind result 



hi central di Mere net result 


l ig. X Overall view ot the computed density contour plots over a spacepiane 
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a) upwind result b) central difference result 

Fig. 9 Spanwise plots of the computed density contours over a spaceplane 
<x/c = \W): M. = 1.5. a = 20° 



j) upwind resuit hi central difference result 

Pig K) Spartwise plots of the computed density contours over a s >aeeplane 


than the central difference scheme. In the 
present upwind scheme where the dissipation 
terms arc constructed mi the matrix torm, each 
characteristic wave has its own minimum dissipa- 
tion. On the other fund, central difference 
scheme where dissipation terms arc constructed 


in the walar torm. requires amount ot dissipation 
which is large enough tor all the waves. This is 
the mam reason that the present upwind is levs 
dissipative and lias better resolution Thus, the 
use ot the proper upwind scheme is recom- 
mended lor vortka! How simulations at high 
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Reynolds number. 

From the solution presented here, another 
important conclusion may be that it is dangerous 
to discuss physics of separated vortical flows 
based on the numerical solutions without suffi- 
cient resolution. It may be necessary to simulate 
the flow field with a sequence of grid cell size 
even with the high-order upwind schemes. 
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